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INTRODUCTION 


Over the last two decades, there has been substantial interest within NASA in the use of lasers for com- 
munication between orbiting spacecraft. In the environment of space, the inherent advantages of lasers (i.e., 
narrow beamwidth and high directivity) are fully exploited and the problem of weather dependent link 
outages vanishes. Prior to their operational deployment in space, however, it is generally agreed that a com- 
prehensive set of experiments must be carried out to verify the performance and feasibility of these systems. 

Currently, the most promising possibility for conducting such experiments centers on the NASA Advanced 
Communication Technology Satellite (ACTS). Experiments are being considered which involve a laser 
system on the ACTS and one or more ground stations; the ACTS-to-ground laser communication link would 
use a near infrared wavelength of about 0.8 micrometer and binary pulse position modulation, with an 
avalanche photodiode detector in the ground station receiver. 

This report investigates the degree to which the communication link may be expected to be degraded by at- 
mospheric turbulence, by calculating the expected increase in bit error probability due to turbulence, in- 
cluding the ameliorating effect of a finite aperture receiver. 

The calculations indicate that, for representative values of the system parameters, the bit error probability for 
a realistic, turbulent atmosphere increases by only a few percent over the value for a hypothetical nontur- 
bulent atmosphere, for a signal strength of less than two hundred incident photons per bit. With the same 
system parameter values, the percent increase in bit error probability becomes greater as the signal strength 
increases, until the increase reaches about thirty percent at a signal strength of four hundred photons per bit. 
For greater signal strengths the increase becomes even larger, but this is not important because the bit error 
probability itself decreases to less than one error in a billion. 

If these results are translated into the increase in signal strength needed to compensate for the effects of tur- 
bulence, it turns out that, to maintain a bit error rate of one error in a million, the signal must be increased 
from the approximately 254 photons per bit needed for this error level for a quiet atmosphere to about 257 
photons per bit for an atmosphere with an average amount of turbulence, a signal increase of about 0.05 db. 
To maintain an error rate of one in ten million, the corresponding figures are about 300 photons per bit for a 
quiet atmosphere and about 303 for an average turbulent atmosphere, a signal increase of about 0.04 db. For 
greater degrees of turbulence the required increase in signal strength will, of course, be greater. For exam- 
ple, if the degree of turbulence, as measured by the variance of the signal fluctuation, is greater than the 
average degree of turbulence by a factor of five, then the signal must be increased by 0.17 db to maintain an 
error rate of one in a million, or by 0.18 db to maintain an error rate of one in ten million. (All these 
figures are for the representative system parameter values referred to above.) 

In the course of the analysis in this report, formulas are developed for the bit error probability as a function 
of signal strength for an avalanche photodiode detector and binary pulse position modulation. Formulas are 
also developed for the mean and variance of the bit error probability, which fluctuates as a result of the 
turbulence-induced fluctuations of the signal strength; because these formulas involve quadratures and are 
time consuming to calculate numerically, approximate formulas are developed which are easily calculated and 
are also sufficiently accurate for system feasibility studies. 
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MATHEMATICAL TREATMENT OF THE EFFECT OF TURBULENCE 
ON THE BIT ERROR PROBABILITY 

Formulas for the Bit Error Probability Statistics 


Fig. 1 below shows diagramatically the space to ground laser communication link being investigated. 
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The transmitter is assumed to be on the ACTS satellite, which is in geosynchronous orbit, but the analysis 
that follows holds for any situation of this type. 
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The signal consists of a series of pulses with binary pulse position modulation. In this modulation scheme 
the data to be sent are digitized, so that each pulse represents a binary digit, or bit (zero or one). For each 
bit to be sent a time interval, called a word, is allocated. Each word is divided into two equal subintervals, 
called slots. If the data pulse is sent in the first slot it represents a zero; if it is sent in the second slot it 
represents a one. If there were no background radiation and no noise in the receiver, there would be no prob- 
lem in sensing which slot contains the pulse. Because of the presence of background radiation and receiver 
noise, plus the fact that the transmitter sends some signal even when it is turned off, each of the slots in a 
given word will contain, at the output of the detector, a non-zero power level. Because of random noise 
fluctuations there will be some appreciable probability, for a weak enough signal pulse, that the slot con- 
taining the signal pulse will contain less total energy than the slot not containing it. Since the detector 
decides whether a zero or a one has been sent by comparing the total energy content of the two slots (by 
means of an integrating circuit, a delay circuit and a comparator), it follows that the detector has some prob- 
ability of making a wrong decision, which may become intolerably large as the signal becomes very weak. It 
is this “bit error probability,” or, in other words, the bit error rate (the ratio of the number of errors to the 
total number of bits sent, for a large number of bits) that will be analyzed in this report, especially the effect 
of atmospheric turbulence on it. 

For the ACTS to ground configuration the beam, which starts at the transmitter in geosynchronous orbit with 
a veiy small diameter, will be approximately 2,000 feet across when it reaches the earth’s atmosphere. The 
receiving aperture (a telescope mirror) will thus intercept only a very small part of the beam; so small that 
the radiant intensity, neglecting turbulence effects, will be essentially constant across the aperture. Tur- 
bulence in the atmosphere will result in intensity fluctuations across the aperture. If the aperture is large 
compared to the distance over which the random fluctuations become uncorrelated, the collection and effec- 
tive summation by the aperture of a number of independent random fluctuations will tend to smooth them 
out; i.e., the variance of the aperture averaged signal will be less than the variance of the signal at a point. 
This is in fact the case in the ACTS to ground configuration. For example, the receiving telescope at the 
Goddard Space Flight Center has a 48-inch diameter, whereas the turbulence fluctuations tend to become in- 
dependent over a distance of a few centimeters. Because a large collecting aperture tends to reduce the size 
of the signal fluctuation, such an aperture averaged signal will result in a smaller increase in the bit error 
probability than would a signal collected by a very small aperture. 

For any particular type of detector, the bit error probability (i.e., the relative frequency of errors) is func- 
tionally related to the signal strength. The signal strength will be defined here as the number of signal 
photons per bit incident on the detector; it will be denoted by s, and will be treated as a continuous variable, 
even though it would appear to be discrete. This is justifiable if the number of photons is regarded as the in- 
tensity of the radiation (which in the classical theory of radiation is a continuous quantity) divided by the 
energy of a photon. It is also mathematically convenient, because the log-normal probability distribution, 
which is generally used to describe the turbulence induced fluctuations, is a continuous distribution. The bit 
error probability E may thus be written as: 


E = f(s;p) (1) 

where the quantity p represents a vector of values of the system parameters. If attention is focused on s and 
no confusion results, p will usually be omitted and f(s;p) will be written simply as f(s), or even E(s), if 
there is no danger of confusing the variable E with the function f. The quantities E and s are here regarded 
as ordinary variables; i.e., E is the bit error probability for a signal pulse of strength s. 
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If the signal strength is regarded as fluctuating randomly (due to atmospheric turbulence), s becomes a ran- 
dom variable with a probability distribution which will be denoted by P(s). The bit error probability also 
becomes a random variable, with a probability distribution which will be denoted by Q(E). The two random 
variables are still related by the function f, so that the probability distribution of E is determined if the pro- 
bability distribution of s is known. The relation between the two distributions is (see, for example, Reference 
1, Chapter V): 


Q(E) = P(s) • |ds/dE| (2) 

provided the function f is monotonic (so that the derivative is nowhere zero). In the present case the function 
f is monotonically decreasing (as will be shown below), so Equation (2) is valid. It is worth noting that the 
time scale of the turbulence fluctuations is measured in milliseconds, so that the turbulence fluctuations con- 
stitute only a very slow drift compared to the data rate, which in this case is measured in hundreds of 
megabits per second. 

Since s = f _1 (E), Equation (2) gives the probability distribution of E in terms of E, in principle. In prac- 
tice, however, the quantities on the right side of the equation may not be explicitly representable in closed 
form. This happens to be the case in this particular problem, so that (2) is of little direct use; even though 
the distribution Q(E) contains complete information about the behavior of the random variable E, this infor- 
mation is not easily accessible. However, the statistical information about E required for this problem is 
mostly contained in the first two moments of the distribution (the mean and variance), and these may be 
represented explicitly in terms of the distribution of s by the following formulas (see any text on probability 
and statistics, such as Reference 1): 


OO 


E = J E P(s) ds 
o 


( 3 ) 


OO 

Var(E) = / (E - E) 2 P(s) ds 


o 


( 4 ) 


The bar over E indicates the mean, and Var denotes the variance. 

The rest of this report will be devoted to developing explicit formulas for the functional relation between E 
and s and for the probability distribution of s, and to developing approximate formulas for the mean and 
variance of E which are easier to calculate than the ones given by the above expressions. 

It is worth noting that the expression for the mean of E given by Equation (3) reduces (as it clearly must) to 
E(s) as the variance of s goes to zero. To see this, one need only note that the integral of 
P(s) is unity (because P is a probability distribution) and that P(s) therefore becomes infinite as the 
variance of s goes to zero. Thus the function P(s), for vanishing Var(s), is a Dirac delta function, so that: 


As Var(s) -*■ 0, 


E - 


; 


E d(s-s) ds = E(s) 


( 5 ) 
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Before proceeding with the development of the formulas given by (3) and (4), alternative forms will be 
developed for the mean and variance of E as series in the moments of the distribution P, in the hope that ap- 
proximate formulas can be obtained by taking only the first few terms of the series. It will turn out that this 
can be done; in fact, quite accurate formulas result from using only the first term of each series. 


Expression of the Formulas as Series of Moments 

If, in Equations (3) and (4) above, the quantities E and (E-E) 2 are expanded as Taylor series about 
s, the series are integrated term by term and note is taken of the facts that 

OO 00 

_f P(s) ds = 1 and J (s-s) P(s) ds = 0 

o o 


(since the integral of a probability distribution is unity and the first moment about the mean vanishes), the 
following series are obtained: 


— 00 e ® (s) 

E = E(s) + I _ Mi 

i=2 i! 


Var(E) = (E(s) — E) 2 + I D (0 (E(5)-E) 2 M 
i=2 i! 


where is the i th moment of P(s) about s: 

OO 

M ; = / (s-s) 1 P(s)ds 


( 6 ) 

( 7 ) 

( 8 ) 


and E (l \ D (l) denote, respectively, the i th derivatives of E(s) and of (E(s)-E) 2 , evaluated at s. 

These series may be somewhat simplified if it is assumed that P(s) is approximately a normal distribution; that this 
is in fact true for this case, since P(s) is a log-normal distribution with a small variance-to-mean ratio, will 
be shown later. With this assumption, P(s) may be written: 

(s-s) 2 

( 9 ) 


P(s) £ 


1 


1 


• e 2 o 2 

\ r 2n~ o s 


It is then clear that all the odd moments Mj vanish. (In fact this will clearly be true for any distribution which is 
symmetric about its mean.) Further, for the normal distribution, the even moments may be represented by the 
general formula (see Reference 2): 


Mi = — 


X 


2 i/2 (i/2)! 


fori = 2,4, 6, ... 


( 10 ) 


Thus, for example, M 2 = <? , M 4 = 3 o s 4 , etc. The series (6) and (7) then become: 
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oo 


( 11 ) 


E = E(5) + I 

i=2,4,6... 


E (i) (s) 
2* (i/2)! 


a 


i 

s 


_ , 00 D<‘> (E(s)-E) 2 

Var(E) = (E(s) — E) 2 + I - i/2 ~ — o> (12) 

i=2,4,6... 2^ (i/2)! * 

For now these series will be left as they are, since no further progress can be made until a specific form is obtained 
for the function E(s). 

Before proceeding to obtain formulas for E(s), it will be shown why the signal fluctuations caused by atmospheric 
turbulence result in an increase in the bit error probability. 


Why Atmospheric Turbulence Increases the Bit Error Probability 

In deriving Equation (6) the function E(s) was expanded in an infinite Taylor series about s. If instead, E(s) is ex- 
pressed as a finite Taylor formula about s, with a second degree error term, it may be written: 

OO 

E = E(s) + 1/2 / E (2) (s + 6 (s) • (s-s)) (s-s) 2 P(s) ds (13) 

0 

where E ® is the second derivative of E(s), and 6 (s) is between zero and one for each value of s, but is otherwise 
completely unknown. From this formula it is clear thatE is greater than E(5) if the second derivative of E(s) is 
positive for all s; in other words: 

A sufficient condition for atmospheric turbulence to cause an increase in the bit error (14) 
rate is that the second derivative of the function E(s) be everywhere positive. 

Intuitively, this can be understood by visualizing the curve of E vs. s, which (as will be shown below) is 
decreasing and highly nonlinear, with a positive second derivative. As the signal fluctuates about its mean, 
the value of E given by E(s) oscillates along the curve, to either side of the point (s, E(s)); but because of 
the fact that the curve is decreasing with s, but more slowly as s increases, the increase in E as the point 
moves to the left is greater in magnitude than the decrease in E as the point moves to the right, so that the 
average value of E is biased upward from its value at s. 


Derivation of Formula for the Bit Error Probability 
Fig. 2 shows a schematic diagram of the detector. 


SIGNAL 

+ *" 
BACKGROUND 


APD 


■0 

4 


T 

/ ( 

o 


) dt 


THERMAL 
NOISE 
CURRENT 
(ZERO MEAN) 


X 

DECISION 


CIRCUIT 


0 OR 1 


INTEGRATOR 


Figure 2. Schematic Diagram of Detector 


6 





A series of signal pulses with some distribution of intensities, which will be assumed to be a log-normal distribu- 
tion arising from atmospheric turbulence, is incident on the avalanche photodiode (APD). (The intensity distribu- 
tion is regarded as already having been partially smoothed by the collecting aperture; this effect will be described 
in the next section.) Since the signal has binary pulse position modulation (BPPM), each signal pulse occupies 
either the first or second slot of a two-slot word. Each pulse consists of the transmitted signal plus background 
noise; also, the slot not containing the signal pulse actually has some signal in it, due to the fact that the transmitter 
cannot be completely turned off during that time slot. 

The incident photons in each pulse give rise to primary photoelectrons which are then randomly amplified by the 
avalanche photodiode; the resulting current, which includes bulk and surface leakage currents, is summed with the 
noise current from an external preamp and is integrated. The integrated outputs for the two slots are com- 
pared and a decision is made as to which of the two slots contains the signal pulse, by choosing the slot 
which gives rise to the greater integrated output. 

The output of the integrator, denoted by x, may be written as follows, for a slot containing a signal pulse: 

T 

x = [(n e + ib T/q)G + i s T/q] + (1/q) / i n (t) dt (15) 

o 


The quantities in this formula are defined as follows: 

n e = number of primary photoelectrons for a slot containing a signal pulse plus background noise; n e is a 
(dimensionless) discrete random variable with a Poisson distribution: 


R(n e ) 


e -IP 


OlP) n « 


n e ! 


(16) 


for n e = 0, 1 ,2,3, . . . ; the mean of the distribution is: 


n e = HP 


(17) 


rj = quantum efficiency of detector (dimensionless) 
p = s -I- b = total number of photons per pulse (signal + background) 

s = number of signal photons per pulse (dimensionless; log-normally distributed continuous random 
variable; fluctuates slowly, on a millisecond time scale, due to atmospheric turbulence) 

b = number of background photons per pulse (dimensionless; assumed constant) 
b =r b T 

r b = rate of background photons (per second; assumed constant) 

T = duration of pulse (seconds; T is the duration of a slot; each word has two slots in BPPM) 

G = gain of avalanche photodiode (dimensionless; discrete random variable with a rather complicated 
distribution) 

i b = APD bulk leakage current (subject to gain; amperes; assumed constant) 
i s = APD surface leakage current (not subject to gain; amperes; assumed constant) 
q = charge of electron (1.602 x 10“ 19 coulombs) 

i n = random noise current from preamp (amperes; continuous normally distributed random variable with zero 
mean) 
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The quantity x is a (dimensionless) random variable which measures the output of the integrator as an amount of 
charge in units of the charge of an electron. Because the random variables ne and G are discrete, while the random 
variable i n is continuous, the probability distribution of x is the convolution of discrete and continuous 
distributions. Adding to this the fact that the gain G of the avalanche photodiode has a rather complicated 
distribution, it would seem that further analysis would be difficult. 

Fortunately, approximations have been developed for this problem which greatly simplify the analysis. It has been 
shown by Gagliardi and Prati (Ref. 3; see also References 4-6, on which Ref. 3 is based) that a random variable 
such as x, which represents the output, integrated over the pulse duration, of an avalanche photodiode with added 
thermal noise, may be replaced with quite good accuracy (under certain restrictions which are generally satisfied 
in practice) by a continuous random variable with a gaussian (normal) distribution. Gagliardi did not explicitly 
take account of the bulk and surface leakage currents of the photodiode in his analysis; Abshire (Ref. 7) has 
modified Gagliardi’s formulas to take account of these effects. The result of these analyses allows replacement of 
the random variable x by the gaussian random variable y; the distribution of y may be written: 

, _ (y~y> 2 

P y (y) = -= e 2 0 2 ( 18) 

y V 2n- Oy 

where the mean y and the variance o “ of y are given by: 

y = [(s + r b T)r, + i b T/q]G + i s T/q (19) 

o 2 y = G 2 [ftG + (1 -ft) (2— 1/G)] [(s + r b T) r, + i b T/q] + i s T/q + (2kT R T)/(q 2 R) (20) 

The first two terms in the formula for the variance arise from the photodiode; the last term is the variance of the 
thermal noise. The quantities in these formulas are defined as follows: 

G = mean gain of the avalanche photodiode (a constant; dimensionless) 

ft = effective ratio of the APD hole and electron ionization coefficients (dimensionless; approx. 0.02 or less 
for silicon devices) 

k = Boltzmann’s constant (1.3806 x 10 -23 joules/K) 

Tr = receiver noise temperature (K) 

R = detector load resistance (ohms) 

Equations (18) - (20) give a good approximation for the probability distribution of the output of the integrator, for 
a time slot containing a signal pulse. The analysis is the same for a time slot which does not contain a signal pulse, 
except that the number of signal photons s is reduced by a factor M called the modulation extinction ratio; M is a 
measure of how well the transmitter succeeds in extinguishing the signal during a no-pulse (zero) time slot. For a 
no-pulse time slot, then, the output of the integrator may be represented by a random variable y 0 , with a proba- 
bility distribution given by Equations (18) - (20), except that the quantity s is replaced by s/M. 

A random variable Y may now be defined which is the difference of the variables y and y 0 : 

Y = y - y 0 (21) 


8 



The probability distribution of the variable Y gives the probability that the output of the integrator for a time slot 
containing a signal pulse is greater than or less than the output for a time slot not containing a signal pulse by some 
amount. (In the absence of background radiation and detector noise, the variable Y would of course have zero 
probability of being negative.) The probability of a bit error is clearly just the probability that Y is negative; i.e. , 
that a time slot with no signal pulse will produce an integrator output greater than a time slot containing a 
signal pulse: 0 

Probability of a bit error = / P y (Y) dY (22) 

— OO 


where P y (Y) is the probability distribution of Y. (It might be noted that the bit error probability, calculated in this 
way, gives exactly the same result as would be obtained from the formula P(l/0) x P(0) + P(0/1) x P(l), where 
P(l) is the probability that a signal pulse is present in a given time slot, P(0) is the probability that no signal is pre- 
sent, P(l/0) is the probability that a signal is observed when none is actually present and P(0/1) is the probability 
that no signal is observed when one actually is present.) The variable Y, being the difference of two gaussian 
variables, is also a gaussian variable. Further, the two variables y and y 0 are independent, because they are 
associated with different time slots. It follows that the mean and variance of Y are given by: 

Y = y - y 0 (23) 


2 2 
°i =°y 


+ o 


y 0 


Thus, the probability distribution of Y is given by the formula: 


P y (Y) = 


V r 2rT (o y 2 + o 2 ) 


V2 


~(Y - (y - y 0 )) 2 

2 (o2 + o^) 


(24) 


(25) 


If one puts into this formula the expressions for the means and variances and forms the bit error probability accor- 
ding to (22), a formula is obtained which can be put in the following form: 


E = 


1 

\f~Tn (c 2 s + c 3 )‘ /2 



- (Y - ci s) 2 


2(c 2 s + c 3 ) 


dY 


(26) 


where, as before, E is the bit error probability, s is the number of photons in a signal pulse, and the quantities c i , 
C2, C3, are constants formed from the system parameters: 


C] = (1 - l/M)nG 


(27) 


c 2 = G 2 \fiG + (1 - p) (2 - 1/G)] (1 + 1/M) rj (28) 

c 3 = G 2 IPG + (1 - p) (2 - 1/G)] [2r b T r\ + 2i b T/q] + 2i s T/q + (4kT R T)/(q 2 R) (29) 


This way of displaying the bit error probability E was chosen to show explicitly the dependence of E on the signal 
strength s, with the system parameters lumped into the three constants. 
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The formula given by (26) may be put into a simpler form by changing variables: 
Y - cj s 


z = 


\/~T (c 2 s + c 3 ) 


Vi 


and noting that 


1 ,°° -z 2 


o 

one thus obtains: 

E(s) = V 2 


J e dz = V 2 ; 


_ a(s) , 2 
1- - J - J e' Z dz 
0 


(30) 


where a(s) = 


Ci s 


\fT (c 2 S + c 3 ) 


>/2 


(31) 


This formula exhibits the bit error probability in a very simple form, involving only the evaluation of an error 
function. 


Some Properties of the Bit Error Probability 

Monotonicity and First Two Derivatives. It is clear from Eqs . (30) and (31) that : 

a(s) increases monotonically with increasing s (32) 

E(s) decreases monotonically with increasing s (33) 

As s varies from 0 to 00 , E decreases from 1/2 to 0 (34) 

The first two derivatives of E(s) are easily computed (note that the definite integral must be differentiated with 
respect to its upper limit); the following formulas are obtained: 


E<'>(s)=- ' 


-a 2 (s) 


'/2 


a (1) (s) 


(35) 


E (2) (s) = 


1V2 


— />2 


a z (s) 2 

• [ a (2) (s) — 2 a(s)*o (1) (s)] 


(36) 
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( 37 ) 


where the superscripts in parentheses denote derivatives; the derivatives of a(s) are given by: 

a (1) (s ) = fi (c 2 s + 2 c 3 ) 

4 (c 2 s + c 3 ) 3/2 

a ( 2 ) (s) = _ c i c 2 (c 2 s + 4 c 3 ) 

8 (c 2 s + c 3 ) 5/2 


(38) 


Note that fi « 1 , while M, G » 1 ; it follows from Eqs. (27) - (29) that Cj , c 2 and c 3 are all positive. From (37) 
and (38) it then follows that a^Ms positive and a®is negative, and therefore from (35) and (36) that E (1) is 
negative and E^is positive: 


(1) (s) < 0 

for all s ^ 0 

(39) 

(2) (s) > 0 

for all s > 0 

(40) 


The fact that the second derivative of E(s) is positive justifies the earlier assertion that atmospheric turbulence 
causes an increase in the bit error rate (see (14) above). 


Dependence on System Parameters . If the expressions for the c’s given by Eqs. (27)-(29) are put into the 
formula for a (s) given by (31), the resulting expression may be cast into the following form: 


a = 


1 _1_ 

s/T (1 - M ) 


” p(5+ i 


-2) + (2— i-> 

G 


1 + 1/M , 2r b T 2i b T 

+ T + 


n s 


n s' 


qq 2 s 2 


II 


'h 


2i s T 


qGVs 2 


+ 


4 k Tr T 


q 2 R G 2 q 2 s 2 


(41) 


From this formula the behavior of a with respect to the various parameters is apparent. Aside from the clear 
monotonic increase of a with s, one sees that a increases monotonically with M, 17 and R, and decreases 
monotonically with fi, r b , T, i b , i s and Tr . If these observations are combined with the fact that E decreases 
monotonically with a, as is evident from Eq. (30), it follows that: 


E increases monotonically with (i, r b , T, i b , i s and Tr ; (42) 

E decreases monotonically with M, rj and R. (43) 


With respect to the mean gain G, the only parameter not yet considered, the behavior of E is more complex. To 
investigate the dependence of a, and thus E, on G, Eq. (41) may be written in the following form: 


a = 


IbG --L 

+ -M 

Vi 

\ G 

G 2 I 



(44) 
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where: a = (l/\/"2) (1 — l/M)>0, sinceM»l; 

c = constant X (1-/3) >0, since fi«l; 
b, d > 0. 

Also, it is clear from Equation (41) that a > 0 for all G and the expression in curly brackets > 0 for all G. 
Equation (44) may be written in the alternate form: 

aG 

{ bG 3 - cG + d } 1/2 

From (45) it follows that: For G = 0, a = 0 
and from (44) that: As G — a -*■ 0 


(45) 

(46) 

(47) 


The first derivative of a with respect to G may be written in the two alternate forms: 

-r [-£-#] 


da 

"dS" 


or 


f b °— §- + w ) 

H-E 


3/2 


G i +cG-2dJ 
{ bG 3 - cG + d ) 


3/2 


From (48) it is clear that: 

As G -*■ °°, da/dG -*■ 0 and d a/dG < 0 

and from (49) it is clear that: 

As G -*• 0, da/dG - a/d' /2 > 0 


(48) 


(49) 


(50) 


(51) 


From (46), (47), (50) and (51), and the fact that a is always positive, it follows that the curve of a vs. G starts at 
zero at G=0, with a positive slope, remains positive, and decreases to zero as but how many maxima 

does the curve have? This question may be answered by setting the derivative equal to zero and noting from (49) 
that this is equivalent to solving the equation: 

f(G) = b G 3 + cG - 2d = 0 (52) 

This is a cubic equation, so it must have either one or three real roots. 
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( 53 ) 


Differentiation of f with respect to G yields: 
df/dG = 3bG 2 + c 

Note that this derivative is everywhere positive, and also note from (52) that f(G) < 0 for G = 0 and f(G) > 0 for 
G-oo ; it follows that the cubic equation has only one real root. Thus a (G) has only one maximum. From the fact 
that E decreases monotonically with a it then follows, from this last result and from (30), (46) and (47), that: 

E -*■ 1/2 as G -*• 0 or 00 


and 


E(G) has only one minimum. 


(54) 


In other words, as G increases from 0 to °°, E decreases from a value of 1/2, goes through a single minimum and 
increases back toward 1/2 (if all the other system parameters and the signal strength are held constant). 


Statistics of the Signal Fluctuation Due to Turbulence, Including Finite Aperture Averaging 


Approximately Log-Normal Distribution of the Aperture- Averaged Signal . In accordance with general prac- 
tice it may be assumed that, at any point of the receiving aperture, the fluctuating signal has a log-normal 
distribution. This may be intuitively understood as resulting from the large number of turbulence-induced in- 
dex of refraction fluctuations over the atmospheric path and the exponential form of the absorption over the 
path: 

L 

- / a (l) dl 

n r = n 0 e (55) 

where n 0 is the number of signal photons per unit area at the top of the atmosphere (the laser beam may be as- 
sumed to be a plane wave, because of the great distance of the transmitter in geosynchronous orbit), n r is the 
number of signal photons per unit area at some point on the receiver aperture (after absorption and scattering over 
the atmospheric path) and a is the total extinction coefficient at the point 1 along the path, where 1 varies from 0 at 
the receiver to L at the top of the atmosphere. Because the extinction coefficient a at each point of the path varies 
randomly due to turbulence, and these random variables are independent except for correlations over a short 
distance, the integral may be regarded as the sum of a large number of essentially independent random variables 
and thus as a random variable itself with a distribution which is approximately normal. The variable n r is 
thus approximately a log-normally distributed random variable. 

The effect of the finite receiving aperture is to sum a number of these (possibly correlated) log-normal variables, 
one for each small area element of the aperture. It is not clear, without further analysis, what the distribution of the 
sum will be. It turns out that the sum (i.e., the total number of photons incident on the aperture, which was called s 
in the preceding analysis, assuming no losses in the receiving optics) is generally approximately log-normally 
distributed (see Refs. 8 and 9). This is surprising, because the central limit theorem states that the distribution of 
the sum of independent random variables approaches a normal distribution as the number of variables becomes 
large, whatever the individual distributions may be. Yet experimental measurements have shown (see Ref. 9) that 
the distribution is log-normal, regardless of the size of the collecting aperture. The reason for this phenomenon is 
discussed in Ref. 9, where it is shown that, because of the skewness and large tail of the log-normal distribution, 
the convergence is slow, so that the log-normal nature of the distribution persists longer than would be expected, 
although of course the distribution will become normal for a large enough number of variables. 
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The distribution of the signal strength s (incident on the detector after aperture averaging) may therefore be 

written in the form: 2 ■> 

_ ('/ 2 In (s/s) + a\T 


p(s) = i 


2 o. 


2 s J~2n s O] 


(56) 


where the quantity o 2 , called the “log-amplitude variance,” is given by: 


o] = 1/4- In (1 + o s 2 /s 2 ) 


(57) 


and s, o 2 are the mean and variance, respectively, of the log-normal distribution of s. In the next subsection 
a formula will be given for the variance o 2 . 


Variance of the Aperture- Averaged Signal . In Ref. 8 a formula is derived (Eq. 6-50, p. 6-24) for the variance 
of the signal intensity after averaging by a finite receiving aperture, assuming plane wave propagation. This 
formula may be adapted to the problem considered in this report by simple changes in units and in some of 
the quantities used; the result is the following formula: 


a 


2 

s 


= 13.1 


s 2 sec 3 9 
A 7/6 


h max 

f (h — h 0 ) 2 C 2 (h)dh 
h 0 


(58) 


where s is the number of photons incident on the detector per signal pulse (considered, as above, as a con- 
tinuous random variable), s is the mean of s, 9 is the zenith angle, A is the area of the receiving aperture in 
square meters, h is the altitude in meters, h 0 is the altitude of the receiver in meters, h max is the altitude (in 
meters) above which the structure constant C 2 is zero and C 2 is the index of refraction structure constant in 
meters _2/3 , assuming the Kolmogorov theory of turbulence. 


The ratio of the standard deviation to the mean is particularly important for the current problem. From Eq. 
58, and recalling that the standard deviation is the square root of the variance, this ratio is given by: 



sec 3/2 9 
A 7/l2 


6 max 

J (h - h 0 ) 2 C 2 

ho 


-i V4 


(h) dh 


(59) 


It is interesting to note that these expressions do not involve the wavelength. This is the case only for large 
apertures; for wavelengths of approximately one micrometer, Eq. (58) will be a good approximation provided 
the aperture diameter is greater than about ten inches (see Ref. 8, page 6-24). The proposed receiver at the 
Goddard Space Flight Center has a diameter of 48 inches; thus the formula should be quite accurate, since 
the proposed wavelength is approximately 0.8 micrometers. 
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Smallness of the Ratio of the Standard Deviation to the Mean . To evaluate Uie integral in Eq. (58) one needs 
to know the dependence of the index of refraction structure constant on altitude. The variation of the struc- 
ture constant with altitude is, in general, extremely complex and it also varies in time. However, a veiy 
simple approximate model for the structure constant is given in Ref. 8 (Eq. 6-20, page 6-14) which is time 
invariant, but (according to Ref. 8) accurate enough to represent a median of the real world. Because of its 
simplicity, this model will be used in some of the following investigations, but there is always the option of 
using more detailed models. The model is defined as follows: 


C» ft) - 



below 20,000 meters above sea level 
above 20,000 meters above sea level 


(60) 


where b = 1.5 x 10“ 13 ; h is the altitude in meters above local ground (assumed to be less than 2500 
meters above sea level) and C n (h) has the same units as given above. (Note that, although h is the altitude 
above local ground, the cutoff of 20,000 meters relates to elevation above sea level.) 

The model of the structure constant may be used to obtain an approximate value for the ratio of the standard 
deviation to the mean, given by Eq. (59). For the receiver at the Goddard Space Flight Center one may use 
the approximate values: 0 = 51 degrees; h 0 = 10 meters; h max = 20,000 meters; A (for the 48" telescope) 

= 1.17 square meters. The integral is easily evaluated to obtain: 


L = 0.03 (61) 

s 

In other words, the standard deviation of the signal strength fluctuation due to turbulence is only about three 
percent of the mean signal strength, for this configuration. The use of a more sophisticated model of the 
structure constant is not likely to qualitatively change this result, so it may be concluded that the fluctuation 
of the signal due to turbulence is quite small compared to the signal itself. 


Normal Approximation to the Log-Normal Distribution . The smallness of the ratio of the standard dev- 
iation to the mean of the fluctuating signal leads to the result that the log-normal distribution given by Eq. 
(56) may be approximated by a normal distribution, with a degree of accuracy which depends on the ratio of 
the standard deviation to the mean. 


To prove this result, one may start by making use of the series expansions: 

1 


In (1 + x) = x — ■ 


x 2 + — L x 2 


-Lx 4 + 

4 


In x = 2 


E x-1 1 / x-1 \ 3 , 1 / x-1 \ 5 , ”1 

x + 1 + 3 \ x + 1 / 5 (x+1 ) ■' J 


for x 2 < 1 ; 


for x > 0. 


(62) 


(63) 
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( 64 ) 


Application of these series to the logarithms appearing in Equations (56) and (57) yields: 
ln(l + y 2 ) = y 2 — j- y 4 +-j- y 6 y 8 + ••• 


(65) 

(66) 


where 


a /s 
s 


In (s/5) = 


= 2 e + £ 3 + £ 5 + ••• 


where 


(s - 5)/s 

£ = 

2 + (s — s)/s 


If one writes (s — s) = no s , where n (not an integer) measures the distance of s from the mean s in units of the 
standard deviation, then: 


2 + ny 

In Eq. (56) there is an s that appears in the denominator; if this s is written in the form s = s + nys and use is made 
of the series expansions above, the log-normal distribution given by Eq. (56) may be written in the form: 


-rr(_ n 2L_) + j_ (_2L_i 3 + ..n +_ld __l 

^ ' 2 + ny / 3 \ 2 + ny / ^ 


D 





P(s) = 



( 68 ) 


If y approaches zero in Eq. (68), a limiting form is obtained by dropping all terms after the first in each series ex- 
pansion, dropping the entire second series in the numerator of the exponential (because it starts with y 2 ) and drop- 
ping the term involving y in the quantity (5 + ny5) in the denominator; the resulting limiting form is; 

- (s - s) 2 

2°s 

W ■ 1 - e (69 

y -*• 0 \f2n o $ 
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which is a normal distribution. It is thus seen that, as the standard deviation becomes small compared to the mean, 
the log-normal distribution approaches the normal distribution with the same mean and variance. However, due to 
the presence of the quantity n in Eq. (68), it is clear that the convergence is not uniform. Because of the smallness 
of the ratio o s /s in this case, as shown in Eq. (61), it may be concluded from the form given by Eq. (68) 
that the log-normal distribution of s may be quite accurately approximated by the normal distribution (69). 


Approximate Formulas for the Bit Error Probability Statistics, Obtained by Truncating 
the Series of Moments 

The fact that the probability distribution of the aperture averaged signal is nearly normal, as shown above, allows 
the use of Eqs. (11) and (12), which were derived on the assumption that the distribution is normal, to develop ap- 
proximate formulas for the mean and variance of the bit error probability E. 

Only the first term of each summation (in Eqs. (1 1) and (12)) will be used; it will be shown below that the resulting 
formulas are in fact good approximations. Taking only the first term of the summation in Eq. (1 1) yields: 

E = E(5) + -j 2) (3) • o s 2 (70) 

If one now takes only the first term of the summation in Eq. (12), substitutes Eq. (70) into the resulting expression 
and simplifies, the result is: 

Var(E) = E (1)2 (5) • a 2 - — E (2)2 (s)-o4 
s 4 s 

It can be shown that the second term of this expression is much smaller than the first; if it is assumed that it is (the 
numerical calculations described in the following section will justify this), the approximate formula for the 
variance may be written: 

Var (E) = E (1)2 (5) . c-. 2 (71) 

The derivatives of E which occur in Eqs. (70) and (71) are given by Eqs. (35) through (38), together with 
Eqs. (27) through (29) and Eq. (31). 

In the next section, a computer program will be described which evaluates the accuracy of the formulas 
given by (70) and (71). 
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COMPUTER PROGRAM TO NUMERICALLY EVALUATE THE ACCURACY OF THE 
APPROXIMATE FORMULAS 


Description of Program 

The computer program, which is written in FORTRAN, evaluates the accuracy of the approximate formulas 
for the mean and variance of the bit error probability E by comparing the values computed from the approx- 
imate formulas against the values computed from the exact expressions given by Eqs. (3) and (4). Listings of 
the program are given in Appendix A. 

The functions E(s) and P(s) occurring in Eqs. (3) and (4) are given by Eqs. (30) (along with Eqs. (31) and 
(27) through (29)) and (56) (along with (57)), respectively. Thus, the exact expressions use the log-normal 
probability distribution, whereas the approximate expressions given by Eqs. (70) and (71) use the normal 
approximation to the log-normal distribution. It should be noted that, although Eqs. (3) and (4) involve no 
approximations, the formula for E(s) given by Eq. (30) is only approximate, because it assumes the approx- 
imate gaussian representation of the photodiode statistics given by Eq. (18); thus even the “exact" values as 
calculated from (3) and (4) are not strictly exact. This should cause no problem, however, because the ap- 
proximation given by (18) is quite accurate. 

The computer program is designed to be used interactively, on a microcomputer. Upon starting the program, 
the user is informed by the program (on the terminal screen) of the current values of the system parameters 
(used in the last run). The user then has the option of changing any or all of these values, along with the 
value of the signal strength. The program then computes the exact and approximate values and prints them 
out, along with the values of the system parameters that were used in the calculation. A sample of the print- 
out is given below in Fig. 3. In this calculation the signal strength was taken as 100 photons per bit, the 
structure constant profile given by Eq. (60) was used and the system parameters used are fairly representative 
values. 

From the results shown in Fig. 3, it is seen that the bit error probability for a non-turbulent atmosphere is 
about 0.0026, that this value increases by about 1.8 percent when turbulence is considered, and that this in- 
crease, as calculated from the approximate formula, is only about 0.2 percent below the increase as 
calculated from the exact formula. For the variance of the bit error probability, the value calculated from the 
approximate formula is about 3.2 percent below the value calculated from the exact formula. The bottom 
three lines in Fig. 3 describe the one standard deviation band about the mean bit error probability; note that 
the bottom of the band (i.e., one s.d. below the mean value of E) lies below the value E(s). (It should be 
noted that, in the computer program, the signal strength s is denoted by the variable N, its mean by NBAR. 
the probability of a bit error by PBE, etc; these variables are clearly defined in the program listings.) 
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EVALUATION OF APPROXIMATE FORMULAS FOR THE MEAN 
AND VARIANCE OF THE BIT ERROR PROBABILITY 


The signal strength is 100.0 photons/bit 

The current values of the system parameters are: 


1. 

Modulation extinction ratio 

= 

100.0 


2. 

Quantum efficiency 

= 

0.7 


3. 

APD mean gain 

= 

250.0 


4. 

APD ionization coeff. ratio 

= 

0.02 


5. 

Background radiation rate 

= 

0.1E+10 

photons/sec 

6. 

Data rate 

= 

0.5E+09 

bits/sec 

7. 

APD bulk leakage current 

= 

0.1E-09 

amperes 

8. 

APD surface leakage current 

= 

0.1E-07 

amperes 

9. 

Receiver noise temperature 

= 

600.0 

degrees K 

10. 

Detector load resistance 

= 

200.0 

ohms 


The current structure constant profile is 1/h. 

The current values of the parameters associated with the signal strength fluctuation are: 


1 . Receiver aperture area 

= 1.167 

sq. meters 

2. Zenith angle 

= 51.0 

degrees 

3. Height of receiver 

= 30.0 

feet 

(above local ground) 



4. Altitude of receiver site 

- 143.0 

feet 

(above sea level) 




Formula for mean 

s^*****************^******************^****^********************* ********** ***************** 


Bit error probability at nbar 
“Exact” increase in BEP 
Approx, increase in BEP 
Error of approx, formula for increase 


0.00257 

1.77880 percent 

1.77471 percent 

— 0.23 percent 


Formula for variance 

**>)<**************************************************************************************** 

“Exact” variance of BEP = 0.241 15E-06 

Error of approx, formula for variance = —3.24808 percent 


Band of one standard deviation about mean BEP 


Mean BEP + 

1 s.d. 

= 0.0031091 

Mean BEP 


= 0.0026180 

Mean BEP - 

1 s.d. 

= 0.0021269 


Fig. 3. Sample Output from Computer Program 
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The computer program uses double precision to calculate the integrals; this is necessary to avoid loss of 
significance due to the sensitivity of some of the calculations. The integrals are computed by subroutine 
DSIMP, which is a double precision subroutine using Simpson’s rule. The user specifies the accuracy re- 
quired; the routine then successively doubles the number of subintervals until the computed value converges 
to the required accuracy. In the calculations an accuracy of one part in a million was generally required, 
with one exception: in computing the table of values of the error function for use by subroutine ERF (which 
computes the error function using either interpolation in the table or a convergent or asymptotic series, 
depending on the value of the argument), an accuracy of one part in ten to the ninth power was required. 

Additional complications in the computation of the integrals in Eqs. (3) and (4) are the infinite limits and the 
narrow structure of the integrand, which latter is due to the small value of the ratio of the standard deviation 
to the mean of E as shown in a previous section of this report. To avoid skipping over the narrow peak in 
the integrand, it was found necessary to require DSIMP to start with a large number of subintervals; a value 
of one million subintervals was used, just to be safe. To handle the infinite limits, the integrals were com- 
puted with successively larger deviations from the mean of the P(s) distribution until convergence was 
achieved; step increases of five standard deviations were used. As might be expected from the preceding 
description, the calculations were quite time consuming; this was not a problem, however, because the 
calculations only had to be done once, to establish the accuracy of the approximate formulas. 

From the results shown in Fig. 3 and from other calculations not shown here, it became clear that the ap- 
proximate formulas are sufficiently accurate for system feasibility studies. In the next subsection the results 
of these calculations are shown as plots of the accuracy of the approximations versus signal strength. A plot 
is also included which shows the increase in the bit error probability due to turbulence, as a function of the 
signal strength (as calculated from the exact formula, using the same system parameter values given in Fig. 
3). 

Computed Results 

Accuracy of the Approximate Formulas . The computer program described above was used to compute the 
accuracy of the approximate formulas for the mean and variance of the bit error probability E, for values of 
the system parameters shown in Fig. 3, which are realistic values. The results are shown in Figs. 4 and 5. 

Fig. 4 shows the percent error of the increase in the bit error probability due to turbulence as computed from 
the approximate formula (70). Note that the quantity plotted is not E, but rather E - E(s), since it is the 
increase that we are interested in, not the value of E itself. From Fig. 4 it is seen that the percent error is 
only a few tenths of one percent for a signal strength of less than 100 photons per bit; for greater signal 
strengths the error increases, reaching about 5.4 percent for a signal strength of 300 photons per bit. Beyond 
this point, although the error continues to increase, the bit error probability itself becomes so small (less than 
10~ 7 ) that the error is not important. 

Fig. 5 shows the percent error of the variance of the bit error probability, as computed from the approximate 
formula (71). The percent error is less than three for a signal strength of under 100 photons per bit, increas- 
ing to about 32 percent at a signal strength of 300 photons per bit. Again, although the error continues to in- 
crease, the bit error probability itself becomes so small as to make the error irrelevant. 


20 



PERCENT ERROR PERCENT ERROR 



AVERAGE SIGNAL STRENGTH IN PHOTONS/BIT 

Figure 4. Accuracy of Approximate Formula for BEP Increase 



Figure 5. Accuracy of Approximate Formula for BEP Variance 
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Increase in Bit Error Probability as a Function of Signal Strength . Table 1 below shows the increase in the 

bit error probability due to turbulence (i.e., E - E(s)), as computed from the exact formula (3), using 
Eq. (30) to compute E(s) and Eq. (56) for P(s); the bit error probability E(s) is also shown. The values of 
the system parameters used are the same as were used above (see Fig. 3). 

Table 1 . Bit Error Probability (BEP) and Increase in BEP due to 
Turbulence as Functions of Signal Strength in Photons 
per Bit, for Representative Values of the 
System Parameters 


Mean Signal Strength s 

Bit Error Prob. 

Percent Increase 

in photons/bit 

E(s) 

in Bit Error Prob 

1. 

0.4754 

0.00014 

2. 

0.4519 

0.00058 

5. 

0.3873 

0.0040 

10. 

0.2980 

0.017 

20. 

0.1751 

0.074 

30. 

0.1026 

0.17 

40. 

0.6016 E— 1 

0.30 

50. 

0.3536 E— 1 

0.46 

60. 

0.2084 E— 1 

0.66 

70. 

0.1231 E— 1 

0.89 

80. 

0.7290 E — 2 

1.2 

90. 

0.4326 E— 2 

1.5 

100. 

0.2572 E— 2 

1.8 

150. 

0.1956 E— 3 

3.9 

200. 

0.1530 E— 4 

6.9 

240. 

0.2020 E — 5 

10.0 

280. 

0.2691 E-6 

13.6 

300. 

0.9846 E-7 

15.7 

360. 

0.4869 E-8 

23.1 

400. 

0.6599 E-9 

29.0 
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Increase in Signal Required to Compensate for the Increased Bit Error Probability . The information in 
Table 1 is plotted in Fig. 6 as the curves labeled “quiet atmosphere” and “turbulent atmosphere (average).” 
The other four curves in Fig. 6 are obtained as follows: division of Eq. (70) by E(s) yields: 

E-E(s) a /J_ Eg)(s,p) \ q2 
E(s) \ 2 E (s,p) / s 

This equation shows that the percent increase in E is approximately proportional to the signal variance, for a 
given value of s (and given values of the system parameters). The last column in Table 1 shows the percent 
increase in E for the signal variance computed from Eq. (58), using the structure constant for an atmosphere 
of average turbulence given by Eq. (60). Because the percent increase in E is approximately proportional to 
the signal variance, as shown by the above equation, the values in the last column of Table 1 may be 
multiplied by any constant to obtain the increase in E for various degrees of turbulence, as measured by the 
signal variance. The top four curves in Fig. 6 were constructed by using degrees of turbulence of two, three, 
four and five times the average amount. 

The increase in the signal needed to maintain the bit error probability at a given level when turbulence is 
taken into account may be read off from the curves in Fig. 6. For example, to maintain E at a level of one 
error in a million (E = 10 -6 ), the signal must be increased from about 254 photons per bit to about 257 
photons per bit, or an increase of about 0.05 db,for an average degree of turbulence. For five times the 
average degree of turbulence, the corresponding signal increase is about 0.17 db. To maintain E at a level of 
one error in ten million (E = 10~ 7 ), the corresponding numbers are about 0.04 db and 0.18 db, respective- 
ly. Thus it is seen that the signal increase required to compensate for the effect of turbulence is not large. 
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BIT ERROR PROBABILITY 



PHOTONS PER BIT 

Figure 6. Increase in Bit Error Probability for Various Degrees of Turbulence 
(For Representative Values of System Parameters) 
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SUMMARY AND CONCLUSIONS 


This report investigates the effect of atmospheric turbulence on the performance of a space to ground laser 
communications link operating in the near infrared, using binary pulse position modulation and an avalanche 
photodiode detector. 

Specifically, the effect of turbulence on the bit error rate is analyzed; general formulas are presented for the 
mean and variance of the bit error rate, which is a random variable, being a function of the signal strength 
which fluctuates randomly due to turbulence. The general formulas are not of much practical use, however, 
because they involve integrals with infinite limits and integrands which contain very narrow peaks, so that 
care must be exercised in numerically integrating them. 

Because of the limited usefulness of the general formulas, approximate formulas were developed which are 
easily computed and sufficiently accurate for system feasibility studies. In the process of developing these 
approximate formulas, a very simple formula was developed for the bit error rate as a function of the signal 
strength, expressed in photons per bit. This formula expresses the bit error rate in terms of a single error 
function, in which the upper limit is a function of the signal strength. 

To assess the accuracy of the approximate formulas, a computer program (written in FORTRAN) was 
developed, which compares the mean and variance of the bit error rate as computed from the approximate 
formulas against the same quantities as computed from the general formulas by numerical integration. The 
results obtained with this program showed that the approximate formulas are very accurate for signal 
strengths of about one hundred photons per bit or less, and reasonably accurate (within about five percent for 
the mean and about thirty percent for the variance) for signal strengths of up to three hundred photons per 
bit. For greater signal strengths the accuracy deteriorates, but the bit error rate itself becomes so small (less 
than one in ten million) that the decreased accuracy, and in fact the increase in the bit error rate, becomes 
unimportant. 

In the course of the calculations to assess the accuracy of the approximate formulas, the bir error rate and its 
increase due to turbulence were also calculated as functions of the signal strength. These results show that 
the bit error rate is a 9 extremely rapidly decreasing function of signal strength, varying from near one half 
for signals of only a few photons per bit to less than one error in a billion for signals of four hundred 
photons per bit (for the representative system parameter values given in Fig. 3). Over the same range of 
signal strength, the increase in the bit error rate due to turbulence varies from less than 0.004 percent to 
about thirty percent. 

Translation of these results into the increase in signal strength needed to compensate for the effect of tur- 
bulence shows that the signal must be increased by about 0.05 db to maintain a bit error rate of one error in 
a million, for an average degree of turbulence, and by about 0.17 db for five times the average degree of 
turbulence (as measured by the signal variance). For an error rate of one error in ten million the figures are 
practically the same. (These numbers apply for representative values of the system parameters.) 

From the foregoing results it may be concluded that, for realistic values of the system parameters and for an 
average degree of turbulence, the increase in the bit error rate due to turbulence is rather mild, not exceeding 
thirty percent for signal strengths small enough to warrant even considering the effects of turbulence. The 
corresponding increase in signal strength required to maintain an error rate of one in a million to one in ten 
million is not more than one or two tenths of a db, even for as much as five times the average degree of 
turbulence. 
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APPENDIX A 

COMPUTER PROGRAM LISTINGS 
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PROGRAM TAPROX 


C 

C This routine tests the approximate formulas for the mean and variance of the bit error probability, for 

C the avalanche photodiode, by checking computed values from the approximate formulas against values 

C computed by numerically integrating the exact formulas. 

C 

C The routine may be run any number of times, by responding “yes” to a prompt. The user may re- 
C specify the system parameters and/or the signal strength before each run. 

C ****************************************************************************************** 

c 


c 


BYTE 

INTEGER *2 
REAL*4 
REAL*8 
EXTERNAL 


REPLY, NARRAY(10,6), BLANK 
ITMLST(IO) 

M, MU, IB, IS, PARAM(IO), NBAR, LOWLIM 
PBARIG, Q, QPREV, PBE, ALPHA, DPBEBR 
PBARIG, PVARIG 


COMMON Cl ,C2,C3, NBAR, SIGSQN, DPBEBR 

COMMON/PRFIL/ HO, IPRTSW 


* 

* 

* 

* 

* 


DATA 


N ARRAY/ 'M','M','G','B','R','B',' I ','I','T','R', 
' \'UV VEYB'/I \'B','S\'RV 


t t 

t f 
9 

t t 


r t <-p t t t t -y t r 

/AY \'R\ 

YA Y 

/Tt / 


ft t t t 

ft ft ft 

9 9 9 

ft ft ft ft 


t f 
9 

ft t t t t 
9 9 

ft ft 

9 9 

ft ft 


t T7 t t 
9 ^ 9 

f/~9t t 

9 ^ 9 


'I 


*************************** *************************************************************** 


Mathematical and physical constants 


Pi 

PI = 3.1415926535 


Boltzmann’s constant 


BK = 1.3806E— 23 ! joules/degree K 

Electron charge 


E = 1.602E-19 ! coulombs 

C 
C 

C Specify the values of the system parameters 

C 

C Modulation extinction ratio (dimensionless) 

C 


M = 100. 
PARAM (1) = M 
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Quantum efficiency (dimensionless) 


MU = 0.7 
PARAM(2) = MU 

APD MEAN gain (dimensionless) 


G = 250. 

PARAM(3) = G 

APD ionization coefficient (dimensionless) 


BETA = 0.02 
PARAM(4) = BETA 

Incident background rate (per second) 


RB = 1.E9 
PARAM(5) = RB 

Data rate (bits per second; one bit per two-slot word) 


BITRAT = 500.E6 

TAU = 0.5 * (1. /BITRAT) 

PAR AM (6) = BITRAT 

APD bulk leakage current (amperes) 


IB = 0.1E-9 
PARAM(7) = IB 

APD surface leakage current (amperes) 


IS = 10. E— 9 
PARAM(8) = IS 

Receiver noise temperature degrees K) 


TREC = 600. 

PARAM(9) = TREC 

Detector load resistance (ohms) 


R = 200. 
PARAM(IO) = R 


Prompt user to change any of these values, if desired 
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100 

110 

* 

* 

* 

* 

* 

He 

He 

He 

He 

He 


111 


c 



122 

123 

1240 
C 

1241 


125 

126 
1242 

C 

127 


C 

C 


TYPE 110, M,MU,G,BETA,RB,BITRAT,IB,IS,TREC,R 
FORMAT(///T2,'The current values of the system parameters are:'// 


T10,'l. Modulation extinction ratio 
T10,'2. Quantum efficiency 
T10/3. APD mean gain 
T10/4. APD ionization coeff. ratio 
T10,'5. Background radiation rate 
T10/6. Data rate 
T10/7. APD bulk leakage current 
T10/8. APD surface leakage current 
T10/9. Receiver noise temperature 


\t43,f8.1/ 

\t47,f7.4/ 

',t43,f8. 1/ 

'447,17.4/ 

' , t46 ,e 1 2 . 4 , t6 1 , ' photons/sec ' / 
',t46,el2.4,t61,'bits/sec'/ 

' ,t46,e 1 2 .4 ,t6 1 , 'amperes 7 
',t46,el2. 4, t61, 'amperes'/ 
',t42,fl0. 2, t61, 'degrees K'/ 
',t42,fl0.2,t61, 'ohms'// 


T9, TO. Detector load resistance 
T2,'Do you want to change any of these values? (Y or N): ',$) 
ACCEPT 111, REPLY 
FORMAT (Al) 

IF (REPLY .EQ. 'N') GO TO 300 


TYPE 121 

FORMAT(//T2,'Type the numbers for the items you wish to change,'/ 
T2,'one at a time, each followed by a carriage return;'/ 
T2, 'follow the last item by a zero and a carriage return:') 
DO 1240 1 = 1,10 

ACCEPT 123, ITMLST(I) 

FORMAT(I3) 

IF (ITMLST(I) .EQ. 0) GO TO 1241 
CONTINUE 


DO 1242 J= 1,10 
JJ = J 

ITEM = ITMLST(JJ) 

IF (ITEM .EQ. 0) GO TO 127 

TYPE 125, (N ARRAY (ITEM , JJ) , JJ = 1 ,6) , PARAM(ITEM) 
FORMAT(T 10, 'Change ',tl7,6Al,t24,' = ',t26,G12.4,T40,'to: ',$) 
ACCEPT 126, P ARAM (ITEM) 

FORMAT(G12.4) 

CONTINUE 

M = PARAM(l) 

MU = PARAM(2) 

G = PARAM(3) 

BETA = PARAM(4) 

RB = PARAM(5) 

BITRAT = PARAM(6) 

IB = PARAM(7) 

IS = PARAM(8) 

TREC = PARAM(9) 

R = PARAM(IO) 

TAU = 0.5 * (1. /BITRAT) 
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Calculate values of the C’s 

300 Cl = (l.-l./M) * MU*G 
C 

BRACKC = BETA*G + (1. -BETA)*(2. - l./G) 

FACT1 = (G**2)*BRACKC 

C2 = FACT1 * (l. + l./M) * MU 

C 

FACT2 = 2 . *RB*TAU*MU + 2.*IB*TAU/E 

TERM2 = 2.*IS*TAU/E 

TERM3 = 4.*BK*TREC*TAU/(R*E**2) 

C3 = FACT1*FACT2 + TERM2 + TERM3 

C 
C 

C Ask user to input value of mean signal strength 

C 

TYPE 301 

301 FORMAT(///T2,'Type mean signal strength in photons/bit (real): ',$) 
ACCEPT 302, NBAR 

302 FORMAT(G12.4) 

C 

C 

C Calculate the variance of the signal strength 

C (sgfluc uses dsimp to evaluate the integral of 

C the stmcture constant over altitude) 

C ^ ^ jjij ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

c 

SIGMAN = SGFLUC (NBAR) 

SIGSQN = SIGMAN* *2 
C 

C Compute PBEbar, the mean bit error probability 

C (use DSIMP; start with an upper limit of mean 

C plus 10 sigma and add 5- sigma’s until integral 

C converges) 

C **************************************** 

c 

UPPLIM = NBAR + 10.*SIGMAN 

LOWLIM - AMAX1 ((NBAR - 10.*SIGMAN), 0.) 

INCR = 10 
Q = 0.D0 
C 

400 UPPLIM = UPPLIM + 5.*SIGMAN 

LOWLIM = AMAX1 ((LOWLIM - 5.*SIGMAN), 0.) 

INCR = INCR + 5 
QPREV = Q 
C 
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402 


c 

c 

c 

c 

c 


403 


* 


ERRBD = l.E-6 
MAXDUB = 14 

CALL DSIMP (DBLE(LOWLIM) , DBLE(UPPLIM), 

ERRBD, 6, MAXDUB, MAXDUB+1, PBARIG, 
Q, INDERR, RELERR, NUMDUB) 


Display intermediate values of the integral for the mean, 
so the user can follow the calculation 


TYPE 403, LOWLIM,UPPLIM,Q, INDERR, ERRBD, RELERR,MAXDUB, NUMDUB 
FORMAT(//T2/PBE mean:'// 

= \T22,F8.2/ 

= ',T22,F8.2/ 

= M24,D22,14// 

= ',T25,I2/ 

= ',T24,E10.2/ 

= '.T24.E10.2/ 

= \T25,I2/ 

= ' ,T25,I2//) 


T14/LOWLIM 
T14/UPPLIM 
T14.TNTGRL 
T14. 'INDERR 
T14, 'ERRBD 
T 14. 'RELERR 
T14. 'MAXDUB 
T14. 'NUMDUB 


IF (INDERR .EQ. 0) GO TO 404 
IF (INDERR .EQ. 2) MAXDUB = MAXDUB + 1 
IF (MAXDUB .GT. 16) GO TO 404 
GO TO 402 

404 IF (INCR .EQ. 15) GO TO 400 
IF (INCR .GE. 30) GO TO 405 

IF (DABS(Q - QPREV)/QPREV .GT. l.D-4) GO TO 400 

405 PBEBAR = SNGL (Q) 

DPBEBR = Q 


Compute the bit error probability for NBAR 
PBENBR = SNGL (PBE(DBLE(NBAR))) 


Compute PBEbar from the approximate formula 

SALPHA = SNGL (ALPHA (DBLE(NBAR))) 

Q1 = SQRT(2.) * Cl 
Q2 = C2*NBAR + 2.*C3 
Q3 = (C2*NBAR + C3)**1.5 
ALPHA 1 = (Q1*Q2)/(4.*Q3) 
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Q2 = Q2 + 2.*C3 

Q3 = Q3 * (C2*NBAR + C3) 

ALPHA2 = ~(Q1*C2*Q2)/(8.*Q3) 

C 

Q4 = ALPHA2 - 2.*SALPHA*ALPHA1**2 
Q5 = EXP( — SALPHA**2) 

PBE2 = — (l.*SQRT(PI)) * Q5 * Q4 
C 

APBEBR = PBENBR + 0.5*PBE2*SIGSQN 


Compute the fractional error of the approximate formula 
for the increase in bit error probability 

IF (PBEBAR .EQ. PBENBR) ERRMN = 0. 

IF (PBEBAR .NE. PBENBR) ERRMN = (APBEBR -PBEBAR)/(PBEBAR- PBENBR) 


Display values of the parameters for this run 


710 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


TYPE 710, NBAR, M,MU,G,BETA,RB,BITRAT,IB,IS,TREC,R 
FORMAT( 

T1 8, 'EVALUATION OF APPROXIMATE FORMULAS FOR THE MEAN'/ 
T18/AND VARIANCE OF THE BIT ERROR PROBABILITY'/ 

T18 ' sjc jjc Hi ifc sjc jf; ^ sjt sj: ifc ^ iji sfs jji ijc jjc jfc }Jc 5jc ijs ijc ^ ^ i}c j$c ij* ^ ijc j|e ^ ^ ^ jje ijt ^ f // 


T2,'The signal strength is ' ,t25,g 12 .4,t39, 'photons/bit'// 
T2,'The current values of the system parameters art:' II 


T10,'l. Modulation extinction ratio 
T10/2. Quantum efficiency 
T10/3. APD mean gain 
T10,'4. APD ionization coeff. ratio 
T10/5. Background radiation rate 
T10,'6. Data rate 
T10,'7. APD bulk leakage current 
T10/8. APD surface leakage current 
T10/9. Receiver noise temperature 
T9, '10. Detector load resistance 


= ' ,t43,f8. 1 / 

= ' ,t47,f7 .4/ 

= ',t43,f8. 1/ 

= ',t47,f7,4/ 

= ' ,t46,el2. 4, t61, 'photons/sec'/ 

= t46,el2. 4, t61 , 'bits/sec'/ 

= ',t46,el2. 4, t61, 'amperes'/ 

= ',t46,el2. 4, t61, 'amperes'/ 

= ', t42,fl0.2,t61, 'degrees K'/ 

= ',t42,fl0.2,t61, 'ohms') 


Display the fractional increase in BEP and the accuracy 
of the approximate formula for the fractional increase 

TYPE 410, PBENBR, 

1 00. *(PBEBAR- PBENBR) / PBENBR, 

1 00. *(APBEBR - PBENBR) / PBENBR, 
100.*ERRMN 
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c 


410 FORMAT(/T2/Formula for mean'/ 

* T2,'**************'// 

* T2,' Bit error probability at nbar 

* t2,' “Exact” increase in BEP 

* 

* t2, ' Approx. " " " 

* 

* t2, ' Error of approx, formula for increase 

* 


= ',t49,gl4.6/ 

= ',t48,gl4.6, 
t64, 'percent'/ 
= ',t48,g,14.6, 

t64, 'percent'/ 
= ',t48,gl4.6, 

t64, 'percent') 


Compute PBEvar, the variance of the bit error probability 
(use DSIMP; start with an upper limit of mean 
plus 10 sigma and add 5— sigma’s until integral 
converges) 

UPPLIM = NBAR + 10.*SIGMAN 
LOWLIM = AMAX1 ((NBAR - 10.*SIGMAN), 0.) 
INCR * 10 
Q = 0.D0 


500 UPPLIM = UPPLIM + 5.*SIGMAN 

LOWLIM = AMAX1 ((LOWLIM - 5.*SIGMAN), 0.) 

INCR = INCR + 5 
QPREV = Q 

ERRBD = l.E-6 
MAXDUB = 14 

502 CALL DSIMP (DBLE(LOWLIM), DBLE(UPPLIM), 

* ERRBD, 6, MAXDUB, MAXDUB + 1, PVARIG, 

* Q, INDERR, RELERR, NUMDUB) 


Display intermediate values of the integral 


C 


C 


TYPE 503, LOWLIM, UPPLIM, Q, INDERR, ERRBD,RELERR,MAXDUB,NUMDUB 
503 FORMAT(//T2,'PBE variance:'// 

',T22,F8.2/ 

',T22,F8.2/ 

', T24,D22.14// 

\T25,I2/ 

',T24,E10.2/ 

',T24,E10.2/ 

',T25,I2 / 

',T25,I2//) 


T 14, 'LOWLIM 
T14, 'UPPLIM 
T14/INTGRL 
T14, 'INDERR 
T 14, 'ERRBD 
T 14, RELERR 
T14, 'MAXDUB 
T14, 'NUMDUB 


IF (INDERR .EQ. 0) GO TO 504 
IF (INDERR .EQ. 2) MAXDUB = MAXDUB+1 
IF (MAXDUB .GT. 16) GO TO 504 
GO TO 502 
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504 


IF (INCR .EQ. 15) GO TO 500 
IF (INCR .GE. 30) GO TO 505 

IF (DABS(Q - QPREV)/QPREV .GT. l.D-4) GO TO 500 
505 PBEVAR = SNGL(Q) 

Compute the variance of the bit error probability 
from the approximate formula 

*******.************************ ********* 

PBE1 = — (l./SQRT(PI)) * Q5 * ALPHA 1 
APBEVR = PBE1**2 * SIGSQN 

Compute the fractional error of the 
approximate formula for the increase 

ERRVR = (APBEVR -PBEVAR)/PBE VAR 

Display the accuracy of the approximate 
formula for the variance and the ratio 
of the standard deviation to the increase 
in the bit error probability 

TYPE 610, PBEVAR, 100.*ERRVR 
610 FORMAT(//T2, 'Formula for variance'/ 

L+mr y *T* *T* *T* *T* ‘T* *T* "T* ’T* *T* *T* *T* ‘T’ 'T' T* / / 

* t2,' “Exact” variance of bit error probability = ',t52,gl4.6/ 

* t2,' Error of approx, formula for variance = ',t52,gl4.6, 

* t68, 'percent') 


Print the results of this run, if desired 
TYPE 650 

650 FORMAT)// 

* T2,'Do you want to print the values 
from this run? (Y or N): ',$) 

ACCEPT 651, REPLY 

651 FORMAT(Al) 

IF (REPLY .EQ. 'N') GO TO 4999 

PRINT 710, NBAR, M,MU,G,BETA,RB,BITRAT,IB,IS,TREC,R 

IPRTSW = 0 
DUMMY = SGFLUC(0.) 

IPRTSW = 1 
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PRINT 410, PBENBR, 

* 100. *(PBEBAR — PBENBR) / PBENBR, 

* 100. *(APBEBR- PBENBR) / PBENBR, 

* 100.*ERRMN 

PRINT 611, PBEVAR, 100.*ERRVR 

61 1 FORM AT(/T2, 'Formula for variance'/ 

* t2, '***************** 7/ 

* t2, '“Exact” variance of bit error probability = ',t52,gl4.6/ 

* t2 , ' Error of approx, formula for variance = ',t52,gl4.6, 

* t68, 'percent') 

C 

PRINT 612, AMIN1 (PBEB AR + SQRT (PBEV AR) , 0.5), 

* PBEBAR, 

* AMAX1 (PBEBAR -SQRT(PBEVAR), 0.) 

612 FORMATS 

* t2,'Band of one standard deviation about mean BEP'/ 

* t2,' **************************************** 7 / 

* t2,'Mean BEP + 1 s.d. = ',t23,gl4.6/ 

* t2,'Mean BEP = ',t23,gl4.6/ 

* t2,'Mean BEP - 1 s.d. = ' ,t23 ,g 1 4.6/////////////) 

C 

C 

C Run another test? 

C ************** 

C 

4999 TYPE 5000 

5000 FORM AT(///T2 , ' Do you want to run another test, for different'/ 

* t2, values of the system parameters and/or the 7 

* t2, 'signal strength? (Y or N) :',$) 

ACCEPT 5001, REPLY 

5001 FORMAT(Al) 

IF (REPLY .EQ. 'Y') GO TO 100 
C 

STOP 

END 

C 

C 

REAL FUNCTION PBE*8 (N) 

c ************************** 

C 

C This routine computes the bit error probability as a function of 

C the signal strength N, for the APD detector. 

c 

REAL*8 N, ERF, ALPHA 

C 

C****************************************************************************************** 

c 

PBE = 0.5D0 * (1.D0 - ERF(ALPHA(N))) 

C 

RETURN 

END 
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REAL FUNCTION ERF*8 (ALPHA) 

*** **** ***************** ******* 

This routine computes the error function, for a given value alpha. 


* 

* 


INTEGER OPENSW 

REAL*8 ROOTPI , TERM , SUM , F ACT 1 , FACT2 , ALPHSQ ,EXPN ,DENOM , 
E(401), E1,E2,E3, ALPHA 1 , ALPH A2 , ALPH A3 , 

TERM 1 , TERM2 , TERM3 , ALPHA 
DATA ROOTPI/1.772453850905516D0/, OPENSW/1/ 


Use different computational methods, 
depending on the value of alpha 
******************************* 

IF (ALPHA .LT. 1.D0) GO TO 1000 

IF (ALPHA .GE. 1.D0 .AND. ALPHA .LE. 5. DO) GO TO 2000 

IF (ALPHA .GT. 5.D0) GO TO 3000 


ALPHA LESS THAN 1; use convergent series 
(first eleven terms) 

*************************************** 


1000 


TERM = ALPHA 
ALPHSQ = ALPHA* *2 
FACT1 = 1.D0 
FACT2 = 0.D0 
SUM = TERM 


DO 1100 1 = 1,10 
FACT1= FACT1 + 2.D0 
FACT2= FACT2 + 1.D0 

TERM = — (TERM * ALPHSQ*(F ACT 1 - 2 . D0))/(FACT 1 *FACT2) 
SUM = SUM + TERM 
1 100 CONTINUE 


ERF = SUM *2 , DO/ROOTPI 
GO TO 5000 


ALPHA between 1 and 5; use quadratic interpolation in a 
table read in from a disk file (pre-computed by DSIMP) 
************************************************ 

2000 IF (OPENSW .EQ. 0) GO TO 2100 
C 


All 



u u 


OPEN (UNIT = 1, NAME = 'VM: ERFTAB.DAT', TYPE = 'OLD\ 

* DISP = 'KEEP', FORM = 'FORMATTED' 

OPENSW = 0 
DO 2010 1=1,401 
READ (1,2011) E(I) 

2010 CONTINUE 

2011 FORM AT(F2 1.14) 

C 

2100 L = IDINT (100. DO* ALPHA) -100 +1 
El = E(L) 

E2 = E(L+ 1) 

E3 = E(L + 2) 

C 

ALPHA 1 = DBLE(FLOAT(IDINT( 100. DO* ALPHA))) * 0.0 IDO 
ALPHA2 = ALPHA1 + 0.01 DO 
ALPHA3 = ALPHA 1 + 0.02D0 
C 

TERM1 = (ALPHA2 - ALPHA)*(ALPHA3 - ALPHA)*E1 
TERM2 = 2 . D0*( ALPHA — ALPHA 1 ) *( ALPH A3 — ALPH A) *E2 
TERM3 = - (ALPHA - ALPHA 1 )*(ALPH A2 - ALPHA)*E3 

SUM = TERM1 + TERM2 + TERM3 
C 

ERF = 0.5D4*SUM 
GO TO 5000 
C 
C 

C ALPHA greater than 5; use an asymptotic series 

C (first eleven terms) 

C 

3000 SUM = 1.D0 

TERM = 1.D0 

FACT1 = — 1.D0 

ALPHSQ = ALPHA**2 
FACT2 = 2. DO * ALPHSQ 
C 

DO 3100 = 1=1,10 

FACT1 = FACT1 + 2. DO 

TERM = - TERM*(FACT 1 /FACT2) 

SUM = SUM + TERM 
3100 CONTINUE 
C 

EXPN = DEXP(- ALPHSQ) 

DENOM = ROOTPI * ALPHA 

ERF = 1.D0 — SUM*(EXPN/DENOM) 


5000 RETURN 
END 
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REAL FUNCTION ALPHA *8 (N) 


This function computes the value of alpha for a given value of 
the signal strength N, for given values of the C’s 

REAL *4 NBAR 

REAL*8 N, ROOT, XNUM,DENOM, ALPHA, DPBEBR 

COMMON C1,C2,C3, NBAR, SIGSQN, DPBEBR 


DC1 = DBLE (Cl) 

DC2 = DBLE(C2) 

DC3 = DBLE(C3) 

ROOT = DSQRT (DC2*N + DC3) 

XNUM = DC1 * N 

DENOM = DSQRT(2.D0) * ROOT 

ALPHA = XNUM/DENOM 

RETURN 

END 


SUBROUTINE DSIMP (XL, XU, ERRBD, BCTMIN, BCTMAX, RACTVT, FCT, 

* SF, ERR IND, ERROR, B COUNT) 

C 

C 

C PURPOSE 

C //////// 

c 

C TO COMPUTE AN APPROXIMATION FOR THE DEFINITE INTEGRAL OF THE 

C FUNCTION FCT(X) (SUPPLIED BY THE USER) BETWEEN THE LIMITS 

C XL AND XU. THE INTEGRATION IS ACCOMPLISHED BY THE REPEATED 

C APPLICATION OF SIMPSON’S RULE, SUCCESSIVELY DOUBLING THE 

C NUMBER OF SUB-INTERVALS UNTIL THE DESIRED ACCURACY IS REACHED. 

C THE SUBROUTINE IS WRITTEN IN DOUBLE PRECISION, TO MINIMIZE 

C ROUND-OFF PROBLEMS. 

C 

C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DESCRIPTION OF ARGUMENTS 

///////////////////////// 


XL - LOWER BOUND OF THE INTERVAL OF INTEGRATION (REAL*8) 

XU - UPPER BOUND OF THE INTERVAL OF INTEGRATION (REAL*8) 

ERRBD - DESIRED RELATIVE ACCURACY OF THE COMPUTED VALUE 
OF THE INTEGRAL (REAL*4) 

(FOR EXAMPLE, IF THE INTEGRAL IS DESIRED TO BE 
ACCURATE TO ONE PART IN A MILLION, SET ERRBD = l.E-6.) 


BCTMIN - THE MINIMUM NUMBER OF SUB-INTERVAL DOUBLINGS 
NECESSARY TO GET A FINE ENOUGH GRID TO HIT ALL 
THE STRUCTURE OF THE INTEGRAND FUNCTION 

BCTMAX - THE MAXIMUM NUMBER OF SUB-INTERVAL DOUBLINGS 
(INTEGER, BETWEEN 3 AND 30, INCLUSIVE) 

RACTVT - A SWITCH TO ACTIVATE THE TEST FOR LIMITATION OF 

ACCURACY DUE TO THE ACCUMULATION OF ROUND-OFF ERROR 
(INTEGER - CAN BE 0 OR ANY POSITIVE INTEGER) 

THE TEST FOR ACCUMULATION OF ROUND-OFF ERROR IS 
ACTIVATED WHEN THE NUMBER OF SUB-INTERVAL DOUBLINGS 
REACHES THE VALUE OF RACTVT. 


THE VALUE OF RACTVT SHOULD BE CHOSEN LARGE ENOUGH SO 
THAT THE APPROXIMATION GIVEN BY SIMPSON’S RULE FOR 
RACTVT +1 DOUBLINGS IS CLEARLY CLOSER TO THE TRUE 
VALUE OF THE INTEGRAL THAN THE APPROXIMATION FOR 
RACTVT DOUBLINGS. (THE PURPOSE HERE IS AVOID 
POSSIBLE SPURIOUS INDICATIONS OF ROUND-OFF ERROR 
ACCUMULATION DUE MERELY TO PECULIARITIES OF THE 
INTEGRAND.) 

IF IT IS NOT EASY TO JUDGE AN APPROPRIATE VALUE FOR 
RACTVT, THE USER MAY SET RACTVT = 0. 

(FOR MOST INTEGRANDS THIS SHOULD NOT CAUSE ANY 
SPURIOUS INDICATIONS.) 


IF THE USER DOES NOT WISH TO APPLY THE TEST FOR 
ROUND-OFF ERROR, HE MAY SET RACTVT TO ANY INTEGER 
GREATER THAN BCTMAX. 


FCT - DUMMY NAME FOR THE DOUBLE-PRECISION FUNCTION SUB- 

PROGRAM USED TO COMPUTE THE INTEGRAND. 

THIS FUNCTION SUB-PROGRAM MUST BE SUPPLIED BY THE 
USER; ITS ACTUAL NAME MUST BE PASSED TO DSIMP BY 
THE CALLING PROGRAM. 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


(NOTE THAT IN THE CALLING PROGRAM THE ACTUAL NAME 
MUST BE SPECIFIED AS THE NAME OF AN EXTERNAL ROUTINE 
BY AN “EXTERNAL” STATEMENT.) 

SF - FINAL APPROXIMATION TO THE INTEGRAL RETURNED TO THE 

CALLING PROGRAM (REAL*8) 

ERRIND - AN INDICATOR WHICH TELLS WHETHER THE ROUTINE DSIMP 
WAS ABLE TO COMPUTE THE INTEGRAL TO THE DESIRED 
ACCURACY (INTEGER) : 

ERRIND = 0 - IT WAS POSSIBLE TO REACH THE RE- 
QUIRED ACCURACY WITHIN BCTMAX 
DOUBLINGS. 

ERRIND = 1 - IT WAS IMPOSSIBLE TO REACH THE 

REQUIRED ACCURACY BECAUSE OF THE 
ACCUMULATION OF ROUNDING ERRORS. 

ERRIND = 2 - THE REQUIRED ACCURACY COULD NOT 

BE REACHED WITHIN BCTMAX DOUBLINGS 
OF THE NUMBER OF SUB-INTERVALS. 
BCTMAX SHOULD BE INCREASED. 

ERROR - THE ESTIMATED MAGNITUDE OF THE RELATIVE ERROR OF SF 
(REAL*4). 

BCOUNT - THE NUMBER OF SUB-INTERVAL DOUBLINGS CARRIED OUT 
DURING THE CALCULATION OF THE INTEGRAL. 

SPECIFICATION STATEMENTS 

//////////////////////////////////////////// 

INTEGER BCTMIN, BCTMAX, ERRIND, BCOUNT, M, S3CT, R3CT, K, 

M2, RACTVT 

REAL*4 ERRBD, ERROR 

REAL*8 XL, XU, FCT, SF, H, SUM, SN, RANGEN, S(3), RANGE(3), 

SEND, SEVEN, SODD 


SF = 0.D0 ! Set integral to zero if the 

IF (XU .EQ. XL) RETURN ! lower and upper limits are 

the same. 
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INITIALIZE 

///////////////// 

BCOUNT = 0 
H = XU-XL 

M = 1 

SEND = FCT(XL) + FCT(XU) 

SEVEN = 0. 

SODD = 0. 

S3CT = 0 
R3CT = 0 

COMPUTE SUCCESSIVE APPROXIMATIONS TO INTEGRAL 

llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllilllllllllllll 

100 BCOUNT = BCOUNT + 1 

H = H/2.D0 

M = 2*M 

S EVEN = S EVEN + S ODD 

SODD = 0. 

M2 = M-l 
DO 101 K= 1,M2,2 

101 SODD = SODD + FCT(XL + K*H) 

SUM = SEND + 2.D0*SEVEN + 4.DO*SODD 

SN = (H/3.D0) * SUM 

IF (BCOUNT .LT. BCTMIN) GO TO 100 

GO TO 200 

1 10 IF (BCOUNT .LT. BCTMAX) GO TO 100 


SF = SN 
ERRIND = 2 

ERROR = RANGEN/DABS(S(3)) 

RETURN 

TEST FOR CONVERGENCE - COMPUTE SUCCESSIVE RANGES 

////////////////////////////////////////////////////////////////////////////////////////// 


200 

C 

c 


c 


c 


S3CT = S3CT + 1 
S(S3CT) = SN 

IF (S3CT .LT. 3) GO TO 110 

RANGEN = DMAX1 (S(1),S(2),S(3)) - DMIN1 (S(1),S(2),S(3)) 
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IF (RANGEN/DABS(S(3)) LT. ERRBD) GO TO 250 
GO TO 300 
C 

210 S(l) = S(2) 

S(2) = S(3) 

S3CT = 2 
GO TO 1 10 
C 

250 SF = S(3) 

ERRIND = 0 

ERROR = RANGEN/DABS(S(3)) 

RETURN 

C 

C 

C TEST GROUPS OF THREE SUCCESSIVE RANGES 

C FOR ROUND-OFF ERROR LIMITATION 

C ///////////////////////////////////////////////////////////////////// 

c 

300 R3CT = R3CT +1 

C 

RANGE(R3CT) = RANGEN 
C 

IF (R3CT .LT. 3) GO TO 210 
C 

IF (RANGE(2) .GE. RANGE(l) .AND. RANGE (3) .GE. RANGE(l) 

* .AND. BCOUNT .GE. RACTVT) GO TO 350 

C 

RANGE(l) = RANGE(2) 

RANGE(2) = RANGE(3) 

R3CT = 2 
GO TO 210 
C 

350 SF = S(l) 

ERRIND = 1 

ERROR = RANGE( 1 )/DABS(S( 1 )) 

RETURN 


END 
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REAL FUNCTION PN*8 (N) 

************************* 

C 

C This routine computes the log-normal probability distribution of the 
C signal strength N, for given values of the mean nbar and variance 
C sigsqn. 

C 

£****************************************************************************************** 


C 

REAL*4 NBAR 

REAL*8 N,Q1 ,SIGSQL,SIGMAL, BRACK, EXPT,DENOM,PI,DPBEBR 

C 

COMMON C1,C2,C3, NBAR, SIGSQN, DPBEBR 

C 

DATA PI/3. 141592653589793D0/ 


********************************************************** **************** **************** 


IF (N .EQ. 0.D0) GO TO 100 

DNBAR = DBLE(NBAR) 
DSGSQN = DBLE(SIGSQN) 


Compute the log-amplitude variance 


C 

C 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Q1 = 1.D0 + DSGSQN/DNBAR**2 
SIGSQL = DLOG(Q 1 )/4 . DO 
SIGMAL = DSQRT(SIGSQL) 


Compute the exponent 


BRACK = 0.5D0*DLOG(N/DNBAR) + SIGSQL 
EXPT = - BRACK**2/ (2 . DO/ SIGSQL) 


Compute the denominator 

DENOM = 2 . D0*DSQRT (2 . D0*PI) *N*SIGM AL 


Compute the distribution 


PN = DEXP(EXPT)/DENOM 
RETURN 
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c 

100 PN = 0.D0 
RETURN 
END 


REAL FUNCTION PBARIG*8 (N) 

This routine computes the product of the functions PBE(N) and 
PN(N), to form the integrand of the integral defining PBEbar; 
this integrand is used by DSIMP and so must be in double 
precision. 

REAL*8 N, PBE, PN 

PBARIG = PBE(N) * PN(N) 

RETURN 

END 


REAL FUNCTION PVARIG*8(N) 

This routine computes the product of the functions (pbe(n)-pbebar)**2 
and pn(n), to form the integrand of the integral defining the variance 
pbevar of the bit error probability; this integrand is used by dsimp 
and so must be in double precision. 


REAL*4 NBAR 

REAL*8 N, PBE, PN, DPBEBR 

COMMON Cl ,C2,C3, NBAR, SIGSQN, DPBEBR 

PVARIG = ((PBE(N)-DPBEBR)**2) * PN(N) 

C 

RETURN 

END 
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REAL FUNCTION SGFLUC*4 (NBAR) 

C 

C This routine uses DSIMP to calculate the integral of the index of 
C refraction structure constant over altitude, to get the variance 
C of the signal strength fluctuation. 

C 

C The user may choose one of several models for the structure constant 

C altitude profile, as well as the values of several parameters. 

C 

INTEGERS PROFIL, PROFL 

C 

REAL *4 NBAR, INTGRL 

REAL*8 PRFL1, XIGRL 

C 

EXTERNAL PRFL1 

C 

COMMON/PRFIL/ HO, IPRTSW 
C 

DATA PROFIL/ 1/, AREA/1.167/, THETA/51./, RCVALT/143./, 

* HO/30./, IPRTSW/ 1/ 

IF (IPRTSW .EQ. 1) GO TO 99 
C 

C PRINTOUT MODE: 

C Print out profile model name and parameter values 

C 

IF (PROFIL .EQ. 1) PRINT 10 

10 FORM AT(/T2, 'The current structure constant profile is 1/h.') 

C 

PRINT 20, AREA, THETA, HO, RCVALT 
20 FORM AT(/T2, 'The current values of the parameters associated'/ 

* t2,'with the signal strength fluctuation are:'// 

* tl0,'l. Receiver aperture area = ', t47,fl0. 3, t65,'sq. meters'/ 

* tl0,'2. Zenith angle = ',t47,fl0. 3, t65, 'degrees'/ 

* tl0,'3. Height of receiver = ',t47,fl0.3,t65,'feet'/ 

* tlO, ' (above local ground)'/ 

* tl0,'4. Altitude of receiver site = \t47,fl0.3,t65,'feet'/ 

* tlO, ' (above sea level)'/ 

C 

SGFLUC = 0. 

RETURN 
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Choose the index of refraction structure constant 
altitude profile and specify values of parameters 

*** ************************************* 


Structure constant profile 


TYPE 100, PROFIL 

FORM AT(///T2, This is subroutine SGFLUC - 7/ 

t2,/The following structure constant profiles are provided: 7/ 
tlO,'l. 1/h profile7 
tl0,'2. ... profile 7/ 

t2,'The current profile is,t25,I2,t27,';7 
t2 , ' if you want to change this, type the7 
t2,' number of the profile desired; if not, type zero: ',$) 

ACCEPT 101, PROFL 

FORMAT(I2) 

IF (PROFL .NE. 0) PROFIL = PROFL 


Area of receiver aperture 


TYPE 110, AREA 

FORM AT(///T2, 'Current area of receiver aperture = ',t39,fl0.4, 
t51,'sq. meters;'/ 
t2,'if you want to change this,'/ 
t2,'type the new value (real);'/ 
t2,'if not, type 0. ; ',$) 

ACCEPT 111, ARA 

FORMAT(G14.4) 

IF (ARA .NE. 0.) AREA = ARA 


Zenith angle of atmospheric path 


TYPE 120, THETA 

FORMAT(///T2, 'Current zenith angle = ',t25,f5.1,t32, 'degrees;'/ 
t2,'if you want to change this,'/ 
t2,'type the new value (real);'/ 
t2,'if not, type 0. : ',$) 

ACCEPT 121, THET 

FORMAT(G10.2) 

IF (THET .NE. 0.) THETA = THET 
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Height of receiver (above local ground) 


TYPE 130, HO 

130 FORM AT(///T2, 'Current height of receiver (above local ground = 

* t55,f8.1,t65,'feet;7 

* t2,'if you want to change this,'/ 

* t2,'type the new value (real)”/ 

* t2, 'if not, type 0. : ',$) 

ACCEPT 131, H00 

131 FORMAT(G14.2) 

IF (H00 .NE. 0.) HO = H00 


Altitude (above sea level) of receiver site (not receiver) 


TYPE 140, RCVALT 

140 FORMAT(///T2, 'Current altitude (above sea level)'/ 

* t2,'of receiver site (not receiver) = ',t27,f8.1,t47,'feet;'/ 

* t2,'if you want to change this,'/ 

* t2,'type the new value (real);'/ 

* t2,'if not, type 0. : ',$) 

ACCEPT 141, RCVAL 

141 FORMAT(G14.2) 

IF (RCVAL .NE. 0.) RCVALT = RCVAL 


Call DSIMP to compute the integral of the profile' 

MAXDUB = 14 

200 IF (PROFIL .EQ. 1) CALL DSIMP (DBLE(H0*0.3048), 

* DBLE(20.E3-RCVALT*0.3048), 

* l.E-6, 8, MAXDUB, MAXDUB+1, 

* PRFL1, 

* XIGRL, INDERR, RELERR, NUMDUB) 

IF (PROFIL .EQ. 2) CALL DSIMP ( ) 

IF (INDERR .EQ. 0) GO TO 250 
MAXDUB = MAXDUB + 1 
GO TO 200 
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Compute standard deviation of signal strength fluctuation 


250 INTGRL 
ROOTIG 
SECTHE 
COEFF 
FACTOR 


SNGL (XIGRL) 

SQRT (INTGRL) 

l./COS (THETA*(3. 14159/180.)) 

3.62 * (SECTHE** 1.5) ' AREA**(7./12.) 
COEFF * ROOTIG 


SGFLUC = FACTOR * NBAR 


Display ratio of standard deviation to mean signal strength 

TYPE 300, FACTOR 
300 FORMAT(/// 

* T2, 'Ratio of standard deviation to mean signal strength = 

* t57,f8.4//) 

RETURN 

END 


REAL FUNCTION PRFL1*8 (H) 


This routine computes the integrand of the altitude integral of the 
index of refraction structure constant, for use by DSIMP 

C**********!|:***Ht!|!******=l<*H<!|!***Ht***!t:****!(:*!|<**!(:!(:s|:*!|:***si<*********s|:******!|:**:i<5|<*****=t::|:*=|:***!|s*!t:=|:s(: 

c 

REAL*4 HO 

REAL*8 H, B, H0MTRS 

C 

COMMON/PRFIL/ HO, IPRTSW 
C 

DATA B/1.5D-13/ 

C 

C ******* *********************************************************************************** 

c 

HOMTRS = DBLE (H0*0.3048) 

PRFL1 = ((H-HOMTRS)**2) * (B/H) 

C 

RETURN 

END 
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PROGRAM ERFTAB 


C 

C This routine computes a table of values of the error function 
C and stores them in the file DY1: ERFTAB.DAT. The Table is computed 
C at intervals of 0.01, from 1.00 to 5.00, to nine significant digits, 

C for use by routine ERF in computing the error function for values 
C between 1.00 and 5.00; series are used for values of less than 1.00 
C or greater than 5.00. 

C ****************************************************************************************** 

c 

INTEGERS EIND, BCOUNT 

REAL*8 ERF(401), X, ERFIGD, Q, ROOTPI 

EXTERNAL ERFIGD 

DATA ROOTPI/1 .772453850905516D0/ 

C 

C ****************************************************************************************** 

c 

C Compute the error function values by using DSIMP 

C (display them as they are computed) 

Q ******************************************* 

c 

X = 0.99D0 
C 

DO 100 1=1,401 
X = X + 0.01D0 

CALL DSIMP (0.D0, X, l.E-9, 15,0, ERFIGD, 

* Q, EIND, ERR, BCOUNT) 

ERF(I) = (2.DO/ROOTPI) * Q 

TYPE 101, X,ERF(I) 

100 CONTINUE 

101 FORMAT(/T15,'ERF(',T19,F10.6,T29,') = \T34,F10.6) 

Open a file in DY1 and write the ERF array to it 

OPEN (UNIT = 1 , NAME = 'DYLERFTAB.DAT', TYPE = 'NEW', 

* DISP='KEEP', FORM = 'FORMATTED', RECORDSIZE=21) 

C 

DO 200 1 = 1,401 
WRITE( 1,201) ERF(I) 

200 CONTINUE 

201 FORMAT(F21 . 14) 

C 

CLOSE(UNIT = 1 , DISP = 'KEEP') 

STOP 

END 
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TABLE OF VALUES OF THE ERROR FUNCTION 


Table of values of the error function, accurate to nine significant digits, precomputed by 
routine DSIMP, for use by routine ERF in computing the error function for values of the 
argument between 1.00 and 5.00, at intervals of 0.01 (series are used for values less than 
1.00 or greater than 5.00). 

The error function is defined by the formula: 


erf (x) = —jr~ f 
V n o 
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1.00 

0.842700793 

1.49 

0.964897865 

1.98 

0.994892000 

1.01 

0.846810496 

1.50 

0.966105146 

1.99 

0.995111413 

1.02 

0.850838018 

1.51 

0.967276748 

2.00 

0.995322265 

1.03 

0.854784211 

1.52 

0.968413497 

2.01 

0.995524849 

1.04 

0.858649947 

1.53 

0.969516209 

2.02 

0.995719451 

1.05 

0.862436106 

1.54 

0.970585690 

2.03 

0.995906348 

1.06 

0.866143587 

1.55 

0.971622733 

2.04 

0.996085810 

1.07 

0.869773297 

1.56 

0.972628122 

2.05 

0.996258096 

1.08 

0.873326158 

1.57 

0.973602627 

2.06 

0.996423462 

1.09 

0.876803102 

1.58 

0.974547009 

2.07 

0.996582153 

1.10 

0.880205070 

1.59 

0.975462016 

2.08 

0.996734409 

1.11 

0.883533012 

1.60 

0.976348383 

2.09 

0.996880461 

1.12 

0.886787890 

1.61 

0.977206837 

2.10 

0.997020533 

1.13 

0.889970670 

1.62 

0.978038088 

2.11 

0.997154845 

1.14 

0.893082328 

1.63 

0.978842840 

2.12 

0.997283607 

1.15 

0.896123843 

1.64 

0.979621780 

2.13 

0.997407023 

1.16 

0,899096203 

1.65 

0.980375585 

2.14 

0.997525293 

1.17 

0.902000399 

1.66 

0.981104921 

2.15 

0.997638607 

1.18 

0.904837427 

1.67 

0.981810442 

2.16 

0.997747152 

1.19 

0.907608286 

1.68 

0.982492787 

2.17 

0.997851108 

1.20 

0.910313978 

1.69 

0.983152587 

2.18 

0.997950649 

1.21 

0.912955508 

1.70 

0.983790459 

2.19 

0.998045943 

1.22 

0.915533881 

1.71 

0.984407008 

2.20 

0.998137154 

1.23 

0.918050104 

1.72 

0.985002827 

2.21 

0.998224438 

1.24 

0.920505184 

1.73 

0.985578500 

2.22 

0.998307948 

1.25 

0.922900128 

1.74 

0.986134595 

2.23 

0.998387832 

1.26 

0.925235942 

1.75 

0.986671671 

2.24 

0.998464231 

1.27 

0.927513629 

1.76 

0.987190275 

2.25 

0.998537283 

1.28 

0.929734193 

1.77 

0.987690942 

2.26 

0.998607121 

1.29 

0.931898633 

1.78 

0.988174196 

2.27 

0.998673872 

1.30 

0.934007945 

1.79 

0.988640549 

2.28 

0.998737661 

1.31 

0.936063123 

1.80 

0.989090502 

2.29 

0.998798606 

1.32 

0.938065155 

1.81 

0.989524545 

2.30 

0.998856823 

1.33 

0.940015026 

1.82 

0.989943156 

2.31 

0.998912423 

1.34 

0.941913715 

1.83 

0.990346805 

2.32 

0.998965513 

1.35 

0.943762196 

1.84 

0.990735948 

2.33 

0.999016195 

1.36 

0.945561437 

1.85 

0.991111030 

2.34 

0.999064570 

1.37 

0.947312398 

1.86 

0.991472488 

2.35 

0.999110733 

1.38 

0.949016035 

1.87 

0.991820748 

2.36 

0.999154777 

1.39 

0.950673296 

1.88 

0.992156223 

2.37 

0.999196790 

1.40 

0.952285120 

1.89 

0.992479318 

2.38 

0.999236858 

1.41 

0.953852439 

1.90 

0.992790429 

2.39 

0.999275064 

1.42 

0.955376179 

1.91 

0.993089940 

2.40 

0.999311486 

1.43 

0.956857253 

1.92 

0.993378225 

2.41 

0.999346202 

1.44 

0.958296570 

1.93 

0.993655650 

2.42 

0.999379283 

1.45 

0.959695026 

1.94 

0.993922571 

2.43 

0.999410802 

1.46 

0.961053510 

1.95 

0.994179334 

2.44 

0.999440826 

1.47 

0.962372900 

1.96 

0.994426275 

2.45 

0.999469420 

1.48 

0.963654065 

1.97 

0.994663725 

2.46 

0.999496646 
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2.47 

0.999522566 

2.96 

0.999971618 

3.45 

0.999998934 

2.48 

0.999547236 

2.97 

0.999973334 

3.46 

0.999999008 

2.49 

0.999570712 

2.98 

0.999974951 

3.47 

0.999999077 

2.50 

0.999593048 

2.99 

0.999976474 

3.48 

0.999999141 

2.51 

0.999614295 

3.00 

0.999977910 

3.49 

0.999999201 

2.52 

0.999634501 

3.01 

0.999979261 

3.50 

0.999999257 

2.53 

0.999653714 

3.02 

0.999980534 

3.51 

0.999999309 

2.54 

0.999671979 

3.03 

0.999981732 

3.52 

0.999999358 

2.55 

0.999689340 

3.04 

0.999982859 

3.53 

0.999999403 

2.56 

0.999705837 

3.05 

0.999983920 

3.54 

0.999999445 

2.57 

0.999721511 

3.06 

0.999984918 

3.55 

0.999999485 

2.58 

0.999736400 

3.07 

0.999985857 

3.56 

0.999999521 

2.59 

0.999750539 

3.08 

0.999986740 

3.57 

0.999999555 

2.60 

0.999763966 

3.09 

0.999987571 

3.58 

0.999999587 

2.61 

0.999776711 

3.10 

0.999988351 

3.59 

0.999999617 

2.62 

0.999788809 

3.11 

0.999989085 

3.60 

0.999999644 

2.63 

0.999800289 

3.12 

0.999989774 

3.61 

0.999999670 

2.64 

0.999811181 

3.13 

0.999990422 

3.62 

0.999999694 

2.65 

0.999821512 

3.14 

0.999991030 

3.63 

0.999999716 

2.66 

0.999831311 

3.15 

0.999991602 

3.64 

0.999999736 

2.67 

0.999840601 

3.16 

0.999992138 

3.65 

0.999999756 

2.68 

0.999849409 

3.17 

0.999992642 

3.66 

0.999999773 

2.69 

0.999857757 

3.18 

0.999993115 

3.67 

0.999999790 

2.70 

0.999865667 

3.19 

0.999993558 

3.68 

0.999999805 

2.71 

0.999873162 

3.20 

0.999993974 

3.69 

0.999999820 

2.72 

0.999880261 

3.21 

0.999994365 

3.70 

0.999999833 

2.73 

0.999886985 

3.22 

0.999994731 

3.71 

0.999999845 

2.74 

0.999893351 

3.23 

0.999995074 

3.72 

0.999999857 

2.75 

0.999899378 

3.24 

0.999995396 

3.73 

0.999999867 

2.76 

0.999905082 

3.25 

0.999995697 

3.74 

0.999999877 

2.77 

0.999910480 

3.26 

0.999995980 

3.75 

0.999999886 

2.78 

0.999915587 

3.27 

0.999996245 

3.76 

0.999999895 

2.79 

0.999920418 

3.28 

0.999996493 

3.77 

0.999999903 

2.80 

0.999924987 

3.29 

0.999996725 

3.78 

0.999999910 

2.81 

0.999929307 

3.30 

0.999996942 

3.79 

0.999999917 

2.82 

0.999933390 

3.31 

0.999997146 

3.80 

0.999999923 

2.83 

0.999937250 

3.32 

0.999997336 

3.81 

0.999999929 

2.84 

0.999940898 

3.33 

0.999997515 

3.82 

0.999999934 

2.85 

0.999944344 

3.34 

0.999997681 

3.83 
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